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ABSTRACT 



A system of two partial differential equations represents the 
transient heat transfer behavior of compact heat exchanger surfaces 
when subjected to a step change in fluid temperature, A solution is 
presented for this system which includes the effects of longitudinal 
thermal heat conduction. Also presented are the solutions for the two 
limiting cases of zero and infinite longitudinal conduction. The 
numerical results have been compared to those of C, P, Howard indicating 
a significant decrease in computational time and an increase in accuracy 
of results. The revised curves of maximum slope of fluid temperature 
versus NTU should be of practical value in the evaluation of heat- 
transfer data obtained by transient testing of compact heat exchanger 
surfaces. An unusual combination of mathematical techniques is presented 
for the solution of a boundary value problem involving partial differen- 
tial equations. The solution combines the application of Laplace trans- 

I 

formation with a numerical technique developed by H. Hurwitz,Jr,, and 
P, F, Sweifel, and adapted by L, A. Schmittroth for the inversion of 
Laplace transforms. This technique greatly expands the number of cases 
to which Laplace transforms may be successfully applied. 
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NOMENCLATURE 



A total heat transfer area, ft . 

2 

A matrix cross sectional area, ft . 

s 

C step increase in fluid temperature, = 

c^ specific heat of fluid, BTU/lb 

c^ specific heat of solid, BTU/lb 

2 

h heat transfer coefficient, BTU/hr ft '^F, 

thermal conductivity of matrix, BTU/hr ft °F* 

Ji total length of fluid flow path, ft, 

dimensionless conduction parameter, = ^s^s^^f^f 

NTU dimensionless heat transfer units, = hA/w_c.. 

» ’ f f 

t dimensionless time variable, = hA0 /W c , 

0 time, hr, 

u solid temperature, °F, 

u^ reference solid temperature, °F, 

V fluid temperature, °F, 

v^ fluid temperature at X = 0, t = 0 -f ; °F, 

v^ reference fluid temperature, °F, 

fluid mass flow rate, Ib/hr, 

W mass of matrix, lb, 

X distance along the fluid flow path, ft. 



X 



dimensionless position variable 



= X . 
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1, Introduction. 



The "single blow" problem refers to the study of the transient heat 
transfer behavior of compact heat exchanger surfaces for the purpose of 
determining their heat transfer characteristics. Interest in compact 
heat exchangers has grown with the development of gas turbine engines. 

A regenerative cycle is desirable to increase the efficiency of the gas 
turbine. Proper design of a compact heat exchanger requires knowledge 
of the heat transfer characteristics of the various materials in numer- 
ous geometric configurations. Until recently it has been necessary to 
design and build the complete heat exchanger in order to determine its 
heat transfer characteristics experimentally. The objectives of the 
single blow studies are to determine the heat transfer characteristics 
from small test sections of the various matrices and to obtain consider- 
able reduction in experimental costs. 

The aim of the experimental testing is to obtain the exit fluid time- 
temperature history subsequent to a step change in the entrance fluid 
temperature. The experimental history must be compared to a theoretical 
time -temperature history to determine NTU, the dimensionless heat trans- 
fer parameter. The heat transfer coefficient, h, or the Colburn J factor 
can then be calculated. 

In 1950 Locke ^9^ outlined five experimental techniques.^ One of the 
simpler of these has become known as the"maximum slope** method. This 
technique compares the maximum slope of the experimental time -temperature 
history to the theoretical maximum slopes for various values of NTU and 
of , the longitudinal conduction parameter. With a cursory inspection, 
one might not appreciate the efficiency of this single point curve matching 

^Numbers in square brackets refer to the bibliography. 



1 



technique. One maximum slope represents a complete theoretical or a 
complete experimental curve. Therefore, it is unnecessary to calculate 
and plot vast quantities of theoretical time- temperature histories. 

This method eliminates the laborious and inaccurate matching of theoreti- 
cal and experimental time -temperature histories, previously required to 
determine NTU. However, the maximum slope technique has two major limita- 
tions, The first problem is the difficulty in analytically solving the 
system of governing differential equations when longitudinal conductivi- 
ties other than zero or infinity are included. The second problem is the 
magnification of experimental error due to curve matching. This problem 
is discussed in detail later. 

Considerable analytical effort has been made to solve various aspects 
of the single blow problem since 1927 when Nusselt first studied it, 
Schumann [l4] developed a solution in 1929 for zero conductivity in the 
direction of flow. His results are for a porous medium such as gravel but 
have been extended to include matrices consisting of metal balls, wire 
screens, or continuous materials for which the effects of longitudinal 
thermal conduction can be ignored. However, as indicated by Mondt and 

Creswick , usually when the matrix is constructed of a continuous 
material in the flow direction, the effects of longitudinal conduction 
must be considered, Creswick outlined a finite difference technique to 
include this effect but his work was not complete enough to cover the area 
encountered in experimental testing. Mondt [[n]] obtained a closed solu- 
tion for the limiting case of — OO , Figure 1 indicates how drasti- 
cally these two cases differ and shows the need for intermediate curves. 
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Maximum Slope 
Figure 1 

In 1948, F. E. Romie, et. al«jj.2^ developed a system of partial 
differential equations governing the transient heat transfer behavior for 
hollow circular cylinders, Creswick arrived at the same set of equations 
for parallel plates. By a finite difference technique, he numerically 
approximated the solution of the system. However, difficulty with con- 
vergence was sometimes encountered. In 1963, C, P, Howard , guided 
by physical considerations, developed an alternate form of the finite- 
difference equations, and he was able to determine the convergence criteria 
and avoid the problem encountered by Creswick, Howard presented a table 
of results and a set of curves of NTU versus Maximum Slope for various 
values of ^ . The table covers a range quite adequate for most experi- 
mental testing. As the values of NTU increased beyond 60, it became very 
difficult to obtain solutions with the finite difference technique due to 
the large computation times involved. 

The objective of the solution presented in this paper is to verify or 
improve the results obtained by C, P. Howard and to avoid large computation 
times as NTU increases. This solution not only avoids the difficulty at 



3 







I 



j 






;l 



large values of NTU, but becomes faster as NTU increases due to the 
corresponding increase in the time at which the maximum slope occurs. 

Referring to the error magnification due to curve matching, C, P. 
Howard showed in Figure 3-A of his Appendix 3, that the error due 
to the experimental determination of maximum slope is magnified by a 
factor of two or greater when the value of NTU is determined by the 
maximum slope technique. The error magnification factor is defined as 
the ratio of percent error in NTU to percent error in experimental 
Maximum Slope. The essential features of the error factor as a function 
of NTU are indicated in Figure 2. 




NTU 

Figure 2 

This figure shows that for the special case, .A = 0, the lowest 
error magnification is between 2.5 and 2 for NTU greater than five. At 
NTU near two there is a peak for which an accurate determination of NTU 
by maximum slope matching is impossible. p{ is non-zero^ the peak 
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shifts to slightly higher values of NTU and the useful range above NTU 
of four decreases as /\ increases • For example, for ^ = .06, an 
error factor of 3, at NTU = 5, is doubled at NTU = 28, Therefore, with 
this technique, one is forced to accept much larger errors as NTU and 

ri increase. To avoid the inherent magnification of error, a modifica- 
tion can be made to the time-temperature curve matching technique with 
the aid of modern computers. When ^ is finite, non-zero it is manda- 
tory to use a computer to calculate the theoretical values of fluid 

2 

temperature, given NTU, and t , It would be elementary to include a 
least square error curve matching technique in the temperature program. 
Then, given three or more experimental values of fluid temperature and 
their corresponding values of t, the computer would compare the calcu- 
lated theoretical temperatures at the given t*s for various NTU*s and 
determine the NTU with the least square error. The anticipated error 
in any particular theoretical value of temperature is + 1% or less. 



**t” denotes the dimensionless time parameter. 
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2 , Mathematical Technique* 

The approach to be presented solved the system of partial differ- 
ential equations for finite, non-zero A by eliminating the variable 
of solid temperature. The result is a third order partial differen- 
tial equation of fluid temperature as a function of time and x (the 
ratio of distance along the matrix to the total length). This equa- 
tion and the necessary boundary conditions were transformed with re- 
spect to time by Laplace transformations to obtain a third order total 
differential equation with respect to x with parameter s* The equation 
was then solved for the transformed fluid temperature. The transformed 
boundary conditions were applied to determine the three arbitrary co- 
efficients which were functions of the parameter s. Because the trans- 
formed solution was very complicated, numerical inversion was used to 
compute the inverse Laplace transform. 

A Gaussian quadrature method, developed by H. Hurwitz, Jr. and P. F. 
Sweifel [[bj for Fourier Transform Integrals, was adapted by L. A. Schmitt 
roth 1^13] for the numerical inversion of Laplace Transforms. This techni 
que avoids the difficulty encountered with alternating, slowly converg- 
ing functions and proved to be essential for the successful inversion of 
this solution. For the two limiting cases = 0 and /) = OO , a direct 
inversion of the Laplace solution was available. 
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3. 



Solution. 



The governing differential equations for the transient heat transfer 

behavior including a finite, non^negat ive, longitudinal conductivity are 
3 

as follows: 



(1) 



(9t 



— (rvr-u) +-2_ 



u 



NTU ax 



( 2 ) 



aor 



Jx = NTU (a-nr) 



The applicable boundary conditions are the following: 



a. 


V(x,0) = 0 


d. au 


.(o,t) = 0 




9X 


b. 


(x,o) = 0 


e, au. 


(l.t) = 0 




9X 


ax 


c. 


v(0,t) = C 







The following is the approach used to solve this set of equations 
for the fluid temperature, v, and the slope of fluid temperature, 
respectively. In a similar manner the solid temperature, u, and other 
desired values can be obtained. 

Solving equation 2 for u results in 



3nr 

a-t 



( 2 a) a = r\r 

NTU 3X 

The partial derivative of u with respect to t and the second partial of 
u with respect to x are determined from 2a, and then substituted into 
equation 1 to give 



( 3 ) 



0 = 



3 2. 

a nr( x.t) 4_ M-rn a/\r(x,t) r£Qj dnrcx.t') 

ax^ a ax 



NTU^ 9r\T(X,t) NTU 

A d± 3 



z 

dnr (x.t) 
ax at 



The equations were developed from a heat balance by Creswick 
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Each term in equation 3 is transformed with respect to t by Laplace 



transformation along with the following boundary conditions; 



a. v(x,0) = 0 

b. -Mr (x,0) = 0 

ax 



The necessary transforms are 

-f C d^nr } - (x, s) 

^ L 9 x^) 9 

OO 



etc. j 

OO /OO sill 

+ j 5 L^J'Cx,t)je dt 



e 




- - or Cx, 0) + 5 [at 



- - f 5 

ax ^ 



Substitution of these values in equation 3 after applying the above bound- 



ary conditions, results in 

_ <5atu,s) , nTU a V ix.s) ntu rs -n^ dnn x .s) nrcx.s) 



The corresponding auxiliary equation is 



(3b) . + iM ru Cs-hl) a — 5—0 

h ^ 



The general solution in the Laplace s-plane for the fluid temperature is 

Ar(x,s3 =: 63(5)6^^^ 



where rl, r2, r3 are the roots of equation (3b). 

The boundary conditions c, d, and e are transformed and then used to 
determine the coefficients Cl, C2, and C3. The transform of C is 

(ir(o,s) --.C . 
s 
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The transform of d Is — — CCjS) — Q , and similarly for e, 

. . 9X 
S) ^ o. 
ax 

From equation 2a the following equation is obtained; 



(5). ,^CX,t) 

ax 



1 aV j- 

NTU dX^ 



3 /r(x,t) 
ax ^ 



Applying boundary condition c to equation 4 gives 



(6). r\rC0,5) - Cl Cs) + GzCs) -t- ^ . 

S 

Rather than apply boundary conditions d and e directly it is convenient to 
use their equivalent form in terms of v« To apply boundary condition d 
it is necessary to transform equation 5 which results in 

(X'jS).. J: — ^ therefore; 

" MTU 3X2- + ^ 

9ll(o,s) _ 1 a^nr (o,s) 9/\r (o.s) _ ^ 

ax NJTU ax^ ax 

When this boundary condition is applied the following equation is obtained 

+/i,C,(:5)+a2.C2Cs)+;73C3(r5)= O. 

MTU 



The application of boundary condition e is similar, resulting in 

( 8 ). J_U?C|e'+n2C2e''^+«.fC3e”0^-^'C'e’'' +A2C2e”‘ +x/ 3 e''’= o, 

MTU 

*2. V 

( 8 a). Define ^ rO 4 ./!.^ j where n = 1, 2, 3. Then the 

^ NTU V 

following matrix equation can be set up using equations 6 , 7 , 8 . 
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1 

O 


— 


i 

R. 


1 1 

Rz. f?3 


0 


C,(5) 

CiCs) 


0 




R,e"’ 


Rzt’'^ Rje”^ 




Cs(s) 



To solve for Cl, C2, C3, Cramer's Rule is applied. The denominator is 



A- [RzR3(e''l-e^' )-f R.R3(e^'-e’'0 + RiR2.(e'’^_e''‘ ^ , 



The three arbitrary coefficients are the following 

1 i. 

O Rs 

0 Rie'"- Rje"* 



C,(s) = 



Ci(sj = 



i C/s 1 

Ri o R-j 



A 



C3(S) = 



t 



i C/5 

R. Ri 0 

R.e"' Rie"' 0 



1 

s 






RAe"-R,R,e 



A 






I LR.Rae"' -R.Rje'^ 



and 



|[R,Ri,e''‘ -R.Rie"'] 



A A 

Substituting Cl, C2, and C3 into equation 4 results in 
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nr(x,sr= — je”''*+R,R3Ce"'-e”0e"^^+ R,Rz(e'“-e'’' )e*^ 

Evaluating equation 9 at x = X results in the following: 



( 10 ) 



ar(l,s) 



) ■hR.i?3(e''‘ - ^ -e'"' ) 



Recall that the Laplace operation Sf(S) — p(+0) corresponds to 

^ F C't') . It follows directly that (1,^:) is equal to 

dt 

the operation S'/^^(1,S) » nr(i^O)— O 

boundary condition <S. . Therefore, 



( 11 ) 






To determine the exit fluid temperature and the slope with respect 
to time of the exit fluid temperature, it is necessary to find the in- 
verse Laplace transforms of equations 9 and 10 respectively. If the 
values of and r^, i = 1, 2, 3; were not functions of s this step would 
be elementary. Unfortunately, the values of R and r are indeed functions 
of s as indicated by equations 8a and 3b respectively. Therefore, for 
this particular case the task of finding an analytic inverse transform 
is nearly hopeless. 

Rather than proceed further on this approach, a computer program to 

calculate the complex roots for given values of s, NTU, and ^ , is 

used. Consequently, this step precludes any possibility of finding an 

analytical solution if one exists. Now, the value of v(l,s) or at 

dt 

(l,s) can be determined for any set of the parameters, C, s, NTU, and , 
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and the Laplace inversion integral may be used to find v(l>t)* Find- 
ing an adequate numerical technique to evaluate this inversion integral 
was no simple matter* 

When Simpson's Rule was applied to evaluate the inversion integral 
using the real (VR) and imaginary (VI) parts of v(l,s), the results were 
grossly in error, the behavior being as described in reference |^6j and . 
The numerical technique described in reference was found to be greatly 
superior to Simpson's Rule for this particular problem* It incorporates 
a technique to accelerate convergence of an alternating series, and uses 
Gauss-Chebyshev formulas for integration* This technique uses either 
the real part (VR) or the imaginary part (VI) and results in two differ- 
ent forms* Since there are no known values for this solution it was 
necessary to use both of these forms to be confident that the correct 
solutions were obtained. 

According to the theory of Laplace inversion, the inversion integral. 



s = G + iy . When this approach was used with Simpson's Rule none of 
the results were independent of G. Recall the form of the inversion 



initially investigated used t = 7 and G of .5, 1 and 2, any error in the 
numerical technique was greatly multiplied giving very poor results. However; 
this indicated that G is an additional parameter that must be determined 
for each new program. Rather than cause difficulty, this parameter proves 
to be essential. For the same G, it was found that one program gave con- 
sistently low results while the other was consistently high. Then by 

4 

Churchill, op, cit,, p, 176* 




, should be independent of G, where 



integral 




. Since the region 
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varying G of both programs the program using VR can be made to agree 
with the results of the program using VI<» This is the approach that 
was used to make the various programs converge to a solution. The 
numerical technique outlined in j^l^ uses Chebyshev's polynomial to fit 
the functior f(s). It was found that this polynomial is very sensi- 

tive to a small change in G. For example, upon changing G from 0.0 to 
0.05, the product (slope t NTU) changed from 1.0877 to 1.1333 an increase 
of 4.19%. Again according to theory, results should be independent of G. 
Therefore, two new programs were written using a Gauss-Legendre quadra- 
ture format to integrate VR or VI. These programs were found to be 
essentially independent of G for reasonable increments of G. For example 
changing G from 0.0 to 0.06, the slope^NTU changed from 1.10096 to 1.09776 
a decrease of 0.29%. Figure 3 indicates a comparison of these two pro- 
grams. On this basis alone, it is felt that the Legendre polynomial is 
more accurate. Unfortunately, in comparison to the programs using the 
Chebyshev polynomial it is considerably slower. Because of this fact, 
the Legendre programs which include an error analysis are being used to 
determine accurate test values for the determination of the best G for 
each program. Once the best G is determined, the speedier Chebyshev 
programs can be used for long data runs. 

The integration technique described in detail in references and 
>3] will be described here as it was applied to this problem. The inver- 
sion integral is the following: 



( 12 ) 



£'X((s)] = m ^ 



/ G "he 60 






ZTTi Jr.: 



G-t-Co 



Then writing f(s) as the sum of the real and imaginary parts results in 
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FIGURE^ 3 






. ^oo 

= zWi j. C^d^) . 



The inversion integral can then be written as 

where the following condition applies: f(s) - f(s). 

It follows directly that 

(12b) ^C6+Ctj) = 

- - ^(0,-itj). 

In the first integral of equation 12a, let y =-yJ and dy = “dy* , then 

(13) f (() - 1^ M 

Applying the relations l2b and 12c, equation 13 reduces to 
On rearranging and simplifying, this becomes 



(12c) 



Gt ^00 

(1^) -f (•<:) - ^i— ^ [[0C&+(-y)cos yt— '-H(fG+cyj5inyfJ dLj 



By definition of the inverse integral, f(-t) = 0, Thus 



On e 



rr 



J n^)(G+-Ly)(:os yi +4^(6fly) Sinyt^ c)(j 

0 L ^ 



A bar over a quantity indicates the complex conjugate. 
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I 



which reduces to 

OO 



j cos Sin Lftc/ij^ 

Two forms of f(t) can be found by expressing the integrals either in 
terms of the real or the imaginary parts. They are respectively 

Gt .OO , 

(16) f(t)= J 0(G^-i<j)CoStjt d<j ^ and 



TT 



(17) 



. I 

! LiJ(G+Lij)sinqt d(j , 






These two equations are the basis of the two programs mentioned earlier. 
In general discussions VR will represent the function (j) cos yt and VI 
the function ^ sin yt. The Gaussian quadrature formula using either 
Chebyshev's or Legendre's polynomial is used to evaluate the integral 
from zero to infinityo Both functions VR and VI are of the general 
form shown in Figure 4 




Figure 4 

If the function were integrated with respect to y and the summation 
stopped at a, there would be significant error. Two operations are 
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included in this solution to reduce the error. The first operation 
is an averaging technique which accelerates convergence. The second 



operation involves integrating each half wave and then summing them. 
To accomplish this we introduce the following change of variable in 
Equation 17. 



Let X 



— ~~ y therefore (j ^ ^ ^ d ^ ^ 



Applying these to equation 17 results in 

(17a) ^Ct)= - ^^/^^[&+^^(X+-L)]slnTT(X+y^)c)x, 

G>t CO -Hoo. 



-4-nrv 



Now let X =. X— ^ , then 



(HM f(t)= - 



M=0 ^ 



^ ( X V /n 4; ^ S I n TT (X V m f ^ ^ 



Upon supressing the prime, we have 



Gt 






(18) f(t)= 2(;.i)4V[G+i|^(^+'«+2)|co-STrxclx 



When equation 18 is written in terms of partial sums, the result is 



(19) 



f (t) = ^ (t) 

— — Z 

IT /yn— ^ GO 7 



where 



(t) - 51 ^ Xm (t) 

= o 



and 
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I 






(19a) 



I.>(t)= COSClTX)cix 



Again referring to > L* A, Schmittroth applied a Gaussian quadrature 
method to evaluate I (t). I 't) can h< approximated as follows: 



j, 1 COS77’(j7 L_ J 



where 



and 



(19c) L( - 2 j -1 

2(2N-M) > 

-uj = ^(- 4j +'"+^) j = +'" + a) 



The number of points used to fit the function is 2 N, Obviously, the 
larger N is, the greater the degree of accuracy one can expect. This 
is the core of the iterative process and an adequate N is desired to en- 
sure reasonable computer time. This point will be discussed in detail 

N 

later. The N coefficients W are the solutions of the system. 



(19d) 



KJ 



2 L 'A/.'^cos 



2M-2 



(2j-l) TT 

2(2N+1) 



- Rp'-?) 

Rp+i) 



where jp = 1, 2, , . • , N, 



5 



Equation 19 uses Chebyshev*s polynomial for the curve fitting. A 
similar form of t) can be developed from Equation 16. As mentioned 
eailier, an integration technique was developed using a Legendre poly- 
nomial for curve fitting. If equation 16 is used to illustrate this 

approach, thr relation ^ is required. Then equation 16 reduces to 

TT 
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( 20 ) f/t) = ['^<^(6+1-^^) Cos 77-X dx. 

-t ^ 

Now, the Gaussian quadrature method can be applied as above* For the 
Gauss -Legendre formulas, a change of variable was applied to equation 17b 
to set up the infinite sum of integrals from zero to one, one to two and 
so on rather than a sum of integrals evaluated from -1/2 to 1/2* 

The Gauss-Chebyshev and the Gauss-Legendre formulas applied to 
the real or imaginary form of the inversion integral produce four differ- 
ent computer programs. Any one of these should give the correct solution, 
but one form may prove to be better for a particular application* 

For zero and infinite longitudnal conduction, the governing differ- 
ential equations reduce to much simpler forms. The solutions for both 
cases are given in Appendix I and are in complete agreement with the 
solutions of Schumman and Mondt for R = 0 and ^ = OO respectively* 
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4. Computer Program. 

The principal objective of the computer program is to solve for the 
maximum slope of fluid temperature at X - * The main program is design- 

ed to determine the slope of the fluid temperature at X which is the 
numerical inversion of equation 11 • 

The subroutine S00TS2 calculates the real and imaginary coefficients 
of the equation 3b given the parameters NTU, , G, and Y. These co- 
efficients are then put into subroutine. T00TS2. 

Subroutine T00TS2 calculates the real and imaginary parts of the 
roots r^^, r^ of equation 3b. It sends its results back to S00TS2 

which arranges them in the proper order and subscripts them appropriately. 
S00TS2 then provides the proper roots as input values to EVAL. The EVAL 
subroutine provides the main program with the values of the real (VR) 

and imaginary part (VI) of ( d-) ^ given the parameters G, 

9t 

IT C, and the roots from S00TS2. 
t 

The main program will call for N of these values (see equation 19c) 
with which it will calculate I (t) (see equation 19). The programmer is 
free to determine the number of I^(t)’s to sum^ where S/yni't)* 

Fifteen to twenty partial sums appear to be sufficient for the function en- 
countered in this problem. Each partial sum calculated is stored and 
serves as an input to subroutine AVER. 

The purpose of subroutine AVER is to accelerate convergence of the 
partial sums by an averaging technique. This technique is discussed in 

reference . The subroutine will take an arbitrary number of partial 

th 

sums and calculate any specified n average, and return it to the main 
program. It also calculates the difference between successive n^^ averages 
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aor 

dt 
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returning it to the main program. The main program then requires this 
difference to be less than an arbitrarily specifiec value (EP). If the 
difference is larger than EP, the main program increases the number of 
partial sums by LD and calls AVER again until the difference is less 
than EP. At this point, the value of the accelerated partial sum is 
accepted and the slof^ is calculated. However, the slope calculated for 
a particular value of dimensionless time is not the value desired. It 
is necessary to search for the maximum slope by varying time. 

Initially, a search using an incremental step of time was used. 

This was eventually abandoned because it would find relative maxima only 
and it was tediously slow if a poor choice were made for a starting time. 
It should be noted that theoretical relative maxima have been found for 
small values of NTU and negative slopes have also been noted. J. M. 
Bannon]^ is concurrently conducting an experimental thesis and has experi- 
mentally found relative maxima and negative slopes for high mass flow or 
small values of MTU. It would be interesting to compare a theoretical 
time temperature curve for his experimental value of with his ex- 

perimental time- temperature curve. 

The purpose of SLOMAX is to find the maximum slope. To avoid the 
problem of choosing a poor starting time, a reasonable range of t is 
chosen. A rule of thumb is that the largest maximum slope will occur at 
a time less than or equal to NTU. SLOMAX uses a parabola curve fitting 
technique with three points. 

The first and last of the three points are program input values of 
time and the middle point is the arithmetic mean. The slopes are calcu- 
lated for the first two values of time and then tested to ensure that the 
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middle slope is larger than the first* This is necessary to obtain the 
maximum rather than a minimum slope* Then the third value of slope is 
calculated and tested to insure that it is less than the middle slope. 

The program has the ability to adjust the arbitrary values of time to 
ensure the middle slope is the largest* Once this basic shape is obtained 
the program calculates the value of time at which the parabola has a maximum. 
This process is repeated using the calculated maximum time and slope as the 
middle point and the two nearest points from the preceeding iteration. 

When the difference of the new value of maximum slope and the previous 
calculated maximum slope is less than EPl, a program input, the slope is 
accepted as the maximum. Then the values of time, slope and input para- 
meters are printed out. 
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TABLE 2» RESULTS OF SLO3 Al^tD SL05 FOR 7 \ = 0.0 
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TABLE 4. RESULTS OF SL02 AND SL04 FUR 7^ = 0»0 
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Discussion of Results, 



Table 1 indicates that programs SL03 and SL05 are not sensitive to 
G and that both programs give the same solution to six decimal places 
for NTU greater than 5, The variation of G from -0*1 to +0.1 is much 
larger than the variations used in the SL02 and SL04 programs which are 
more sensitive to G, All results presented were calculated with eight 
point formulas* 

One way to increase the accuracy of the program is to increase the 

number of terms in the partial sum. There is a test in the program 

which will increment the number of partial sums if the nth average has 

-9 

not converged to within 1 x 10 of the preceeding nth average for each 
calculation of slope. The test for the sufficient number of partial 
sums indicates that for NTU greater than five, twenty partial sums will 
give the same accuracy as 100 partial sums for programs SL03 and SL05. 

At NTU of five, the SL03 program automatically increased the number of 
partial sums to 25 which then gave as accurate a solution as 95 partial 
sums. There is an indication that for lower values of NTU, more partial 
sums may be necessary. 

R. E, Maximjl^ presented a table of values of maximum slopes cal- 
culated with Schumann® s solution for the special case of = 0,0. 
Schumann’s solution is an infinite series of Bessel functions. Each 
value of slope was calculated using Bessel functions accurate to ten 
significant places and the summation was continued until there was no 
change in the tenth decimal place of the solution. Unfortunately, the 
search by Maxim could not distinguish between relative maximum slope 
values and his search ended when the slope at the next increment of time 
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was less than the preceeding slope o Relative Tnaximum slopes have been 
found by the theoretical solution presented in this paper and their ex- 



Figure 26, The results of Table 2 clearly indicate that Maxim's values, 
though accurate to six places, were calculated at a time much too low to 
have found the largest relative maximum slope in the NTU range less than 
five. For values of NTU equal to and greater than five, both the time 
at which maximum slope occurred and the values of maximum slope times NTU 
calculated by SL03 agree to at least seven decimal places with Maxim's 
results* SL05 solutions were in complete agreement with the results of 
SL03o 

The computation time depends on the degree of accuracy required to 

calculate one value of slope* The accuracy required of the search routine 

-9 

for the results of Table 2 was + 1 x 10 difference in the new value of 
maximum slope as compared to the previous iteration* For values of NTU 
less than ten, the number of iterations required was four or five times 
the number for NTU values greater than ten* Once the approximate value 
of time at maximum slope is known, the search routine can be set to find 
the maximum with fewer iterations* But, even without optimum starting 
times the computation time was approximately half a minute for the 7\ = 0*0 
case* For non zero values of » the computation time per search is 

expected to be twelve minutes* 

Results of SL05 for ^ -*04 are shown to be insensitive to G in 
Table 3* One test point at NTU - 10 indicates good agreement between 
SL05 and SL03 programs* As for the ^ = 0 case, results indicate good 
agreement with C* P* Howard's results * For NTU greater than five 

Howard's results are consistently low but all have less than 1% differ- 
ence as compared to SL05. As in the p{ 0 case, Howard's results for 



istence is supported by experimental evidence found by J* M* Bannon 
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NTU equal to two are much lower than those obtained by SL05. Again this 
is probably a failure of his search technique to distinguish between 
relative maximum slopes. Unfortunately, he did not publish the times 
at which the maximum slopes occurred, Howard also used the incremental 
step search and the search ended when the newly calculated slope was 
less than the preceeding slope. 

Programs SL02 and SL04 are quite sensitive to parameter G, The 
results obtained graphically and listed in Table 4 were compared to the 
accurate solution of Maxim for f{ =0.0, The percent error of the graphi- 
cal results was less than + 1.5% except at NTU = 2, Figure 6 indicates 
that this intercept was not determined. Considerable effort would be 
required to find this intercept and obviously the mean value at NTU = 2 
has unacceptably large error. For the other values of NTU, however, 
intercepts were easily found as indicated by Figure 5, 

Table 5 shows similar results for ^ = 0.005. No accurately known 

values are available for this region, C. P. Howard's results are com- 
pared to the results obtained by SL02 and SL04. If as expected, the 
SL02 and SL04 results are 1 to 1.5% too large, then C. P, Howard’s 
solutions for = 0,005 would be approximately ,5% too low for NTU 

greater than 10, Figure 7 again shows how the intercept or best G was 
obtained. 
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7. 



Conclusions and Recommendations* 



Program SL05, using the VR, Gauss-Legendre formulas, is a fast 
accurate tool capable of supplying all the theoretical maximum slope 
data needed for experimental testing of compact heat exchanger surfaces* 
The SL03, VI, program is equally accurate but slightly slower* The SL05 
program is being applied to verify the NTU versus Maximum Slope curves 
presented by C* P* Howard [sj • The results of these calculations will 
be published at a later date. 

The combination of Laplace transforms with numerical inversion should 
be applicable to most boundary value problems in which one is able to 
express the solution as an inverse Laplace transform* Therefore, this 
technique is primarily limited in application to linear systems of part- 
ial differential equations with constant coefficients* 

The two main problems encountered in the numerical inversion were 
the slow convergence of the integration and radical behavior of the func- 
tion of s at time equal zero* The VR and VI functions were cyclic in 
nature with slowly decreasing magnitudes and period* For this particular 
solution the Legendre polynomial approximated the functions better than 
Chebyshev*s polynomial. However, for other functions a number of other 
polynomials may give better results* The acceleration of convergence is 
adequately accomplished by the present numerical inversion but it gener- 
ates a problem in addition* Because of the averaging technique applied 
to the partial sums to accelerate convergence the formal error analysis 
normally applied to each partial sum is impossible* Therefore, the only 
way remaining to determine the accuracy of the calculated maximum slopes 
was by comparison with R, E. Maxim's accurately determined values for the 



34 



limiting case of ^ = 0. Even though direct inversion of the Laplace 
transform was possible for this case, the numerical inversion was ac- 
complished in seconds with the same accuracy attained by the approxi- 
mation to an infinite sum of Bessel functions carried out to six decimal 
place accuracyo For non-zero 7) , there were no accurate solutions with 
which to compare. For these values an indication of accuracy was obtain- 
ed by calculating the same points with an independent program. The SL03 
program is not independent of SL05 in a strict interpretation of indepen- 
dence. However, they are different programs, SL03 using the imaginary 
part of the inversion integral and SL05 using the real part. Either 
program should give the same result. The same result to six or seven 
decimal places was obtained with these programs indicating excellent 
accuracy for a numerical process. The comparison of SL05 results with 
Howard’s results, obtained with a rather crude finite-difference tech- 
nique, indicated that he obtained very accurate results by judicious ap- 
plication of the finite-difference technique. Most of his results were 
less than one half of one percent low when compared to SL05, 

The following recommendations indicate avenues of investigation 
that could be undertaken with either computer programs developed in con- 
junction with this thesis or the mathematical technique presented. 

It is recommended that the temperature programs TEMP4 and TEMP2 be 
applied to obtain a t ime- temperature history for comparison with corres- 
ponding experimental results to determine how well the mathematical model 
represents the behavior of the experimental test rig. 

An investigation should be conducted to determine the cause of and 
physical significance of the relative maximum slopes obtained both experi- 
mentally and theoretically for the heating cycle. One objective would be 
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to determine if the relative maximum obtained experimentally correspond 
in magnitude and time of occurance with those obtained by theory. Ex- 
perimentally it was noted that the relative maximums are practically 
indistinguishable for the cooling of the solid, but that for the same 
matrix and mass flow of fluid the relative maximums are very distinct 
for the case of heating the solid. An investigation should be made to 
determine theoretically the effect of longitudinal conduction for the 
cooling cycle as compared to the present solution for the heating cycle. 

The objective of this investigation would be to determine if the theoreti- 
cal solution is capable of predicting a difference in response on cool- 
ing as compared to heating. The cooling cycle solution would require 
appropriate sign changes in the heat balance from which governing partial 
differential equations were derived. This may, therefore, require a new 
solution for cooling since the sign changes could cause a change in the 
auxiliary cubic equation 3b of Section 3. It is possible that a sign 
change of a term in this equation could cause a change in the significance 
of the parameter p\ on the result of fluid temperature or slope of fluid 
temperature. Bannon [^l^ has noted experimentally that the largest maximum 
slope obtained by cooling agreed fairly well with that obtained by heating. 

The theoretical results, in the region of NTU = 2, are difficult to 
obtain, and the determination of NTU given an experimental maximum slope 
in this NTU region is impossible. An attempt should be made to develop 
a curve matching technique employing a least square error method. This 
method would use three or more values of time and their corresponding 
temperatures from an experimental time-temperature history. The TEMP4 
program would calculate theoretical values of temperature at the given 
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times for an estimated value of NTU, Then an NTU search routine could 
calculate the sum of the squared error of each temperature « After re- 
peating this operation for two bracketing values of NTU, a Lagrange 
iteration method would calculate the NTU with the least squared error* 
Another possible means of curve matching would be to match the dimension- 
less time at which the maximum slope occurs for the region of NTU = 2, 
rather than matching the magnitudes of maximum slope* This would require 
an adequate starting mark on the experimental time-temperature history 
and a correction of the time at the maximum slope to account for the 
system time delay of the experimental response* 

A more general system of governing partial differential equations 
can be obtained from the heat balance as derived b> Creswlck, hy in- 
cluding a term to represent the energy stored in the fluid* This term 
is required for transient heat transfer testing using a liquid rather 
than a gas. The result of this addition changes equation 2 of section 
3 to the following 



9nr _ 1 

3X MTU 






, the dimensionless fluid capacitance, 
3 



where ~ ^ ^ ^ 

Cs 

and 'y is the density of the fluid In Ib/ft* This equation combined 
with Equation 1 of Section 3 would then be the new governing partial 
differential equations* The present mathematical technique should be 
capable of solving this system of equations* 
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APPENDIX I 



Solution for Zero and Infinite Longitudinal Conduction. 

1. ZERO LONGITUDINAL CONDUCTIVITY 

As in the solution for finite conductivity the development of this 
case is the same initally; cf, Equation 3a, page 8. 

(1) c^AT(x,s) h 

a <5)X a 

For this case where p\ is defined as a, it is necessary to multiply by a. 
The Substitution of a = 0 results in the equation 



(2) r\r' -f b s__ AT = o, 
s+i 

The general solution of equation 2 is 



AT(x,S) 



C,(s) = 



C,(s) 



- fc> 5 V 

e 



The boundary conditions are the same as in the previous solution. 
Boundary conditions ciand b were applied in order to arrive at equation 
1, above. Boundary condition C is v(0,t) = C which transforms into v(0,s) 
C/s. Applying boundary condition C to equation 3 results in 
Ar(O^S)=— = CiCs) ^ therefore 

>s 



ic) 5 

w ^Cx,s)- nr(±>s) = C.e~T+i 



It follows directly, 
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(5) and 3)t 



/7>^-hA 

r C e 5+-I 



This special case was solved analytically in 1929 by T, E, W. 
Schumann • The following development will show that the present 

solution agrees with Schumann's, 



Then 



Let s* = s + 1 and apply it to equation 4. 

^ X 

AT(X, s'-i) -_C_ e s- 

5 -i. 



nr(x^s-l) =, ± 



CO 

w-i s'''^ 



From page 229 of reference [2^ , we find 



or 



(7) 






Comparing equation 6 and 7, the inverse transform of equation 6 becomes 



From page 294 [2j is 

jC"^[f(S-a)] = [f (S)] therefore, 

(9) ru-(x,t)= (i)~ , 

(9a) nrcx,t) - 
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TJ. •• • 




Define a = btx and W = 2 , then compare equation 9a with the form 



(9b) 






Note that 









=■ t 



d 2- 



Equation 9a then becomes 



nr(;x,t-) = (’e 






/>7 - 0 



Schumann’s solution is 



(11) m-CM.a') e"^'^£ g'" d"' (ZiY^) 

J\n — ^ / t t •rt \ 



/V) =:0 



c/(y2.) 



Equations 10 and 11 are equivalent which can be shown as follows, 
page 392 modified Bessel function identify is 





d fvV 




] = V/'' 








Then, for 


n = 0 and ’ 


W = 2 






H 

L_J 

, 'x. 






Apply the 


chain rule 


such that 






dCY] 


- dr 


dw 




d a 


d w ‘ 


d <x 


It follows 


that 








(J [x„ 




II 




da 





where 



From 



fY\- O 
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By mathematical induction it can be shown that 



( 12 ) =a-^I,r,(2{s:). 
da'^ 

Equation 10 then becomes 

(13) tV^[lj2Va)J 
From page 392 ^4j , 

j = L’'^Jm('2i Vo?)^ 

= JUziVST). 

Note that Schumann's development is for the case of heating the fluid as 
time increases and that the present case is developed for cooling of the 
fluid with increasing time, which changes the sign of t in the exponential 
term It then follows that equation 13 is, indeed, 

( 11 ) , 

ck<)£r 

Since the program for numerical inversion of equations 5 and 6 is already 
available it is much easier and probably much faster to apply it than to 
evaluate equation 11. The same theoretical solution was obtained by G. L. 
Locke in 1950 [^9j . 
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2. INFINITE LONGITUDINAL CONDUCTIVITY 



When ^ = CO » the solid temperature is no longer a function of x. 

^ ^ 1 I { I 

In fact the following is true: — ) — 0* Therefore, equation 

o X 

1, page 7, becomes indeterminate and must be replaced by:^ 



I [Arrx,t)-U(0] dx 

( 2 ) ^^0^/0 _ NTU ['u.(t)-nr(x,t)l 



and equation 2 is 



The boundary conditions to be applied to this case are: 
a. v(x,0) = 0, 
be v(x,0) = 0, and 
c. v(0, s) = C/s. 

Solve equation 2 for u(t) and take the derivative of u with respect 
to t. Then substituting these back into equation 1 and applying boundary 
condition c, results in 



di_ ^ ^ Tx.t) j 9Ar Cy,t) — 



MTU 



at 



r/u-(i,t)-ci . 



_ -1 



K/“ri I 



Multiply through by NTU and rearrange the terras to arrive at 



3 r\rC X)i) MTU aar(x,t) . - C- O. 

3xat 9t 



Transformation of each term in 3 by Laplace gives 
^For further discussion, see page 47. 
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(^) ^ru-Cx.o) S^nr(x.S) _|_NTUr~nr^X,o) 4-Snr(x,5) -|-/v(i,5)- 
~ dy 9X ^ -* 

Applying boundary conditions a and b and dividing by s, equation 4 reduces 
to 



(5) i^T-y ar(x.5) = ^ _ r\r(i,s) 

3 ^ 5 



Again by treating s as a parameter, equation 5 is a total differential 
equation. The solution of the differential equation is 



Ar(x,3) X Cx(s) 



^ - /OrCi; 
-S 



Now apply boundary condition c to equation 6 which results In 

0> = I = C, fs) + [f- - 



( 8 ) 



-1 


1- ^ 


s 1 


NTU 




NTU *X 


e~ 


-f - 



+ 



1 



and 

ISTU -5 y 



Solve equation 8 for v(l,s' as follows: 
It A 



NTU 









fi ^1 


_ e 




) -NTU 

- +e 


NTU . S 


s 


HTU . S 
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Define 



5 



; then ^ 



= 1 - j“ NTiJ 



, and 



/ -tJTd) 

(i-e J 



Nru 



- MTU 



(9) 



at(i,s) = S. 

s 



a ?A) 



The particular fluid temperature of interest is that temperature at x = 1* 
Also the value of C must be dimensionless, for comparison with llondt's 
solution, and can be set equal to one. It follows directly that 



ar(i,s) = 1 (i 

(s+ f) 



ATtl,S) = 



J. f-NTU 

5-hf S-h^ 



-h 



f 

s(si-V 



» 



The inverse Laplace transform of equation 10 is 

-NTU‘^e -tl-e" 



which reduces to 
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I 




I 

I 



where 



(11) nr(i.t) = MTU Je 



-ft 



/^ru 



. The solution as a function of 



t ime i s 



nr{±,it) = ^ 7717 — ^ ^ 



This solution is in complete agreement with that derived by J. R. 
Mondt Qll] in 1961. Since a convenient inverse transform was available 
it was unnecessary to apply the intricate numerical inversion. 

To derive Equation (1) for this special case it was necessary to 
refer to the heat balance from which Creswick derived the governing 

equations. The last term of Creswick* s Equation (1) was eliminated and 
the remaining terms were arranged in terms of dimensionless parameters. 
Then both sides of the equation were integrated over the length of the 
matrix to arrive at Equation (1) on page 44. 
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APPENDIX II 



Flow diagrams and program listings. 

The flow diagrams will preceed the applicable program listings. 
The symbols used in the flow diagrams are defined as follows; 



c 






Terminal 



Input/Output 




O 



a 



Decision 



Connector 



Offpage Connector 



the beginning, end, or point of 
interuption in the program 



input symbols are preceeded by 
a terminal. 



branches are labeled by sign 
according to the sign of the 
value enclosed. 

enclosed numbers refer to the 
prefix numbers of the corres- 
ponding statements on the pro- 
gram listing. 

designates an entry to or exit 
from a page. 



There are six distinct programs which are used separately, although 
they could be combined into a larger package to make use of one or another 
at the user's option. There are the four "slope” programs, SL02, SL03, SL04 
and SL05 which have been described in detail in the body of the thesis, and 
the two "temperature” programs TEMP2 and TEMP4, which have been referred to 
in the text but not described in detail. The latter are intended to pro- 
vide a time-temperature history 5 no use of these programs has been made in 
this thesis. See page 35 for recommendations concerning their use. 
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PR0GRA14 SL03* 
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PROGRAK SL03 

D I KENS I ON VR(2C) ,VI (20) ,SUI 200) ,AV( 1 2 , 1 0 ) , R ( 1 29 ) , X { 1 29) ,ER(200) , 
1RR( 126) ,RI ( 123 ) ,T{ 10) , SLO( 10) 

COKMCN RR.RI ,T,SLO,G, R, A,PI ,SL0PE,')T,NT1 ,NT2,EP1 ,11 ,T1 
COKKCN NT3,TK 

101 READ 100,G,B,A,T1,TM,EP1,£P,L1,LD,I 1 
100 F0RMAT{7F10.5,3I3) 

C SET 11=0 ON LAST DATA CARD , 

C II IS THE AVERAGE AND TFE POLYNOMIAL ORDER 

IF( I 1 ) 10,10, 102 

102 CONTINUE 
NT = -1 

. NT2=0 
K=1 

NT1=-1 

PI=3. 1U15926536 
WRITE OUTPUT TAPE i*,301 

301 F0RMAT(24H PROGRAM IS SLOPE 1,'N=8^/) 

WRITE OUTPUT TAPE 4,302 

302 F0RKATI30H EVALUATION SUBROUTINE 2 SHORT/) 

WRITE OUTPUT TAPE 4,305 

305 F0RMAT(25H USES LEGENDRE POLYNOMIAL/) 

204 L=L1 
NT3=-1 
SU=0. 

D=0.0 

E=1.0 

DO 203 J=1 ,100 
CALL GAUSSC(D,£,VAL) 

D=D+1 .0 
E=E+1.0 

SU(J)=VAL+SU( J-1 ) 

503 IF( J-L)203,206,206 
206 L=L+LD 

CALL AVER(SU,AV,I1,J,DIF,SUM) 

500 IF(0IF-EP)202,202,203 
203 CONTINUE 

CALL AVER(SU,AV, 8 , 1 00, Cl F , SUM ) 

202 EG=EXPF(G*T1 ) 

EG2=2.*EG/T1 
SL0PE=EG2*SUK 
PRINT 3,G,B,A,T1 
PRINT 4,SUM,DIF,J,SU( J) 

PRINT 2, Tl, SLOPE 

221 CALL SLOMAX 
TMl=8.*TM/3. 

IF(NT2)406.406, 101 

406 IF(T1-TM1)204,204,222 

222 PRINT 3,G,R,A,T1 

2 F0RMAT(3H T= F 1 0 . 5 , 3X , 6HSL0P£=E20 . 8/ ) 

3 F0RMAT(3H G= Fl 0 . 5 , 3X , 2H 6=F 10. 5, 3X,2 H A=F 1 0. 5 , 3X , 2HT=F 1 0. 5/ ) 

4 FCRMAT(5H SUP= E20. 8 , 3X , 4H0 IF = F 1 0. 9, 3X , 3HSU ( I 3, 2H ) =E20 .3 / ) 

GO TO 101 

10 STOP 10 
END 
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PI 

II 
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li 



101 

100 



102 



PROGRAM SL05 

DIMENSION VR(20),VI(20) ,SU(200),AV( 1 2 , 1 0 ) , R ( 1 29 ) , X ( 1 29) 
1RR( 128) ,RI (128) ,T{ 10) fSLOI 10) 

COMMON RR.RI » J » SLO ,G t B , A , P I , SLOPE T , NT 1 ,NT2, Epl 1 1 1 , T 1 
COMMON NT3,TM 

READ 100,G,B,A,Tl,TM,EPl,EPfLl,LD,I 1 
F0RMAT(7F10.5,3I3) < 

SET 11=0 ON LAST DATA CARD 

II IS THE AVERAGE AND THE POLYNOMIAL ORDER 

IF( II no, 10, 102 

CONTINUE 



,ER(200) 



NT = -1 
NT2=0 
K=1 

NT1=-1 

PI=3. 1415926536 
WRITE OUTPUT TAPE 4,301 

301 F0RMAT(24H PROGRAM IS SLOPE 5, N=8/ ) 

WRITE OUTPUT TAPE 4,302 

302 F0RMATI30H EVALUATION SUBROUTINE 2 SHORT/) 

WRITE OUTPUT TAPE 4,305 

305 F0RMAT(25H USES LEGENDRE POLYNOMIAL/) 

204 L=L1 
NT3=-1 
SU=0. 

D=0.0 
E=1 .0 

DO 203 J=l,100 
CALL GAUSSt)( D,E,VAL) 

D=D+1.0 

E=E+1.0 

SU(J)=VAL+SU( J-1 ) 

503 IF(J-L)203,206,206 
206 L=L+LD 

CALL AVERISU.AV, II .J,DIF,SUM) 

500 IF (DIF-EP)202,202,203 
203 CONTINUE 

CALL AVERISU.AV, 0 , 1 00, D IF , SUM ) 

202 EG=EXPF(G*T1 ) 

EG2=2.*EG/T1 
SLCPE=EG2*SUM 
PRINT 3,G,B,A,T1 
PRINT 4,SUM,DIF,J,SU( J) 

. PRINT 2,T1 ,SLOPE 

221 CALL SLOMAX 
TMl=8.*TM/3. 

IF(NT2)406.406. 101 

406 IF(T1-TM1)204,204,222 

222 PRINT 3,G,B,A.T1 

2 F0RMAT(3H T=F 10.5, 3X, 6HSL0PE=E20.8/ ) 

3 F0RMAT(3H G=F10 .5, 3X, 2HE = F 10 . 5. 3X,2 HA=F 10. 5 , 3X , 2HT=F 1 0. 5/ ) 

4 F0RMAT(5H SUM=E20.8 , 3X , 4H0 I F=F 1 0.9, 3X , 3HSU( 13, 2H ) =E20.8 / ) 
GO TO 101 

10 STOP 10 
END 



f 



/ 



) 
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PROGRAI^i SL02’' 



Input ^ 




5 G, B, A, T1,’TM 



Set II = 0 
On last 
data card 
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PROGRAM SLC2 

DIMENSION RR( 128) ,RI ( 1 28 ) , V I ( 20 ) , ( 20 ) ,SU(200) , AVI 1 2 , 1 0 ) , R ( 1 29 ) , 

1X( 129) ,SLO( 10) ,T( 10) 

COMMON SL0PE,SL0,T,NT,NT1,NT2,K,EP1 tGtB,A,PI,T1,TM 
COMMON NT3 

101 READ 100,G,B,A,T1,TM,TD,EP,L1,LD» II 
100 FCRMAT(7F10.5,3I3) 

C SET II = 0 ON LAST DATA CARD 

IF( 1 1 ) 10, 10, 102 

102 C 1 = 1.0 

210 WRITE OUTPUT TAPE 4,301 

301 F0RMATI24H PROGRAM IS SLOPE 2, N=4) 

WRITE OUTPUT TAPE 4,304 

304 FORMAT! 26H USES CHEBYSHEV POLYNOMIAL) 

401 WRITE OUTPUT TAPE 4,302 

302 F0RMATI30H EVALUATION SUBROUTINE 2S /) 

W1C=. 10776070/. 98480775 

W2C=. 083333333/ .86602540 
W3C=. 045908434/. 64278761 . 

W4C=. 01 2997531/. 342020 14 

NT = -1 

NT1=-1 

NT2=0 

K=1 

EP1 = . 000001 
NT3=-1 

204 P=3. 1415926536/Tl 
DO 500 1=1,200 
. 500 SU(I)=0.0 
L=L1 

Pl=1./9. 

SGN=-1 . 

SU = 0. 

DO 203 J=1 ,100 
603 DO 200 1=1,8 
Y=P1»P 

CALL S00TS2 ( G, B, A, Y , R , X , RR , R I ) 

CALL EVAL2S! RR,RI,G,B,Y,C1 ,VR,VI) 

VR(I)=VR 

VI(I)=VI 

200 P1=P1+(l./9.) 

605 SU( J)=SGN*(W1C«(VI (4)+Vl!5) 1+W2C*!/ I( 3)+VI( 6) )+W3C*( VI! 2)+VI(7) )+ 
1W4C*(VI( 1 )+VI(8) ) )+SU( J-1 ) 

IF(J-l) 212,212,213 

212 PRINT 3,G,B,A,T1 

213 CONTINUE 
SGN=-SGN 

Pl = Pl+( 2./9. ) 

IF ( J-L)203,206,206 
206 L=L+LD 

CALL AVER(SU,AV,I1 ,J,DIF,SUM) 

IF(DIF-EP)202,202,203 
203 CONTINUE 

CALL AVER(SU.AV,10,100,DIF,SUM) 

202 EG=EXPF(G»T1 ) 

EG2=2.*EG/T1 

SL0PE=EG2*SUM 

PRINT 4,SUM, DIF, J,SU( J) 

PRINT 2, T1, SLOPE 
'221 CALL SLOMAX 
T>ll=8.*TM/3. 

IF(NT2)406,406, 101 
406 IF(T1-TM1)204,204,222 
222 PRINT 3,G,B,A,T1 

2 F0RMAT(3H T=F 1 0 .5 , 3X , 6HSL0PE = E 1 5. 8/ / ) 

3 F0RMATI3H G = F 1 0 . 5, 3X, 2HB=F 1 0 . 5, 3X ,2 HA=F 1 0. 5 , 3X . 2HT=F 1 0. 5/ ) 

4 F0RMAT(5H SUM*E 1 5 . 8, 3X , 4HD I F*F1 0 . 9, 3X , 3HSU (13, 2h ) =E 1 5.8 / ) 

GO TO 101 

10 STOP 10 
END 
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PROGRAM SL04 

V n I MENS I ON RR( 128) .RK 128) ,VI (20) , (20) ,SU(200) , AV( 12,10),R(129), 
1X( 129) ,SLO( 10) ,T( 10) 

COMMON SLOPE, SLO,T, NT, NT1,NT2,K,EP1 ,G,B, A, PI, T1,TM 
COMMON NTS 

101 READ 100,G, B,A, T1 ,TM, TD,EP,L1 ,LD, II 
100 FORMAT(7F10.5,3I3) 

C SET 11=0 ON LAST DATA CARD 

IF( 1 1 ) 10, 10, 102 

. 102 Cl = 1 .0 . 

210 WRITE OUTPUT TAPE 4,301 

301 FORMAT! 24H PROGRAM IS SLOPE 4, N=4) ! 

WRITE OUTPUT TAPE 4,304 

304 F0RMAT(26H USES CHEBYSHEV POLYNOMIAL) 

401 WRITE OUTPUT TAPE 4,302 . - 

302 F0RMATI30H EVALUATION SUBROUTINE 2S / ) 

W1C=. 10776070/. 98480775 

W2C=. 083333333/. 86602540 » 

W3C=. 045908434/, 64278761 ■ ' 

W4C=. 01 2997531/. 3420201 4 

NT = -1 

NT1=-1 

NT2=0 

K=1 

EP1=. 000001 
NT3=-1 

204 P=3.1415926536/T1 
DO 500 1=1,200 
500 SU(I)=0.0 
L=L1 

Pl=1./9. 

SGN=1. 

SU=0. 

DO 203 J=1 ,100 
IF(J-1 1600,600,601 

600 P1=-7./18. 

GO TO 603 

601 IF( J-2)602,602,603 

602 Pl=11./18. 

603 DO 200 1=1,8 
Y=P1*P 

CALL S00TS2 ( G, B, A , Y , R, X, RR , R I ) 

CALL EVAL2S(RR,RI,G,B,Y,C1,VR,VI) 

VR( I )=VR 
VI ( I )=VI 

200 Pl=Pl + ( 1 ./9. ) 

IFU-l 1604,604,605 

604 SU( 1 )=.5*SGN*(W1C*(VR{4 )+VR(5) )+W2C*( VR( 3)+VR(6) )+W3C*( VR(2) +VR (7) 
1 )+W4C*(VR{ 1 )+VR(8) ) )+SU( J-1 ) 

GO TO 212 

605 SU( J)=SGN«(W1C*{VR(4)+VR(5) ) +W2C* (VR ( 3 ) +VR ( 6 ) ) +W3C» ( VR( 2) + VR(7) ) + 
1W4C»(VR( 1 )+VR(8) ) )+SU( J-1 ) 

IF(J-l) 212,212,213 

212 PRINT 3,G,B,A,T1 

213 CONTINUE 
SGN=-SGN 

Pl = P1 + (2./9. ) 

IF( J-L)203,206,206 
206 L=L+LD 

CALL AVER(SU,AV, II ,J,OIF,SUM) 

IF(DIF-EP)202,202,203 
203 CONTINUE 

CALL AVER! SU,AV, 10, 100, OIF, SUM) 

202 EG=EXPF(G*T1 ) 

EG2=2.*£G/T1 

SL0PE=EG2«SUM 

PRINT 4,SUM,DIF,J,SU( J) 

PRINT 2,T1 , SLOPE 

221 CALL SLOMAX 
TMl=8.*TM/3. 

IF(NT2)406,406, 101 

406 IF(T1-TM1 1204,204,222 

222 PRINT 3,G,B,A,T1 

2 F0RMAT{3H T=F 10. 5, 3X , 6HSL0PE=E1 5. 8/ / ) 

3 F0RMATI3H G=F1 0 . 5, 3X, 2H8=F 1 0. 5. 3X ,2 HA=F 10. 5 , 3X , 2HT=F 1 0. 5/ ) 

4 F0RMAT(5H SUM=E 1 5 . 8, 3X, 4H0 1 F=F 1 0. 9, 3X , 3HSU ( 13, 2H ) =E 1 5.3 / ) 

GO TO 101 

10 STOP 10 
END 
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PROGRAM TEMPU* 

DIMENSION RR(128),RI(12e),VI(20),VfU20),SU(2Q0),AV( 1 2 ,1 0 ) , R ( 1 29 ) 
1X( 129) ,ER(200) ^ 

COMMON GtB,A,T,TM,TD,EP,Ll ,LO, n,Ml ,RR,RI,PI, 

101 READ 100,G,B,A,T,TM,TD,EP,L1,LD, I 1,M1 
100 F0RMAT(7F10.5,3I3, I 1 ) 

SET 11=0 ON LAST DATA CARD 
IF( 1 1 ) 10, 10, 102 

102 CONTINUE 

PI=3. 1415926536 
WRITE OUTPUT TAPE 4,300 
300 FORMAT! 30H PROGRAM IS TEMPERATURE 4, N=8/) 

WRITE OUTPUT TAPE 4,302 

302 F0RMAT(30H EVALUATION SUBROUTINE 1 SHORT/) 

WRITE OUTPUT TAPE 4,305 
305 F0RMATI25H USES LEGENDRE POLYNOMIAL/) 

204 L=L1 
L=L1 
SU=0. 

0 = 0.0 

E=1.0 

DO 203 J=l,100 

CALL GAUSSQ(C,E,VAL,ER) 

D=D+1.0 

E=E+1.0 

SU( J)=VAL+SU(J-1 ) 

IF! J-L)203,206,206 
206 L=L+LO 

CALL AVER!SU.AV^ II ,J,DIF,SUM) 

IF!DIF-EP)202,202,203 
203 CONTINUE 

CALL AVER!SU,AV, 8, 1 00, C I F , SUM ) 

202 EG=EXPF!G*T) 

EG2=2.*EG/T 
TEMP= EG2*SUM 

WRITE OUTPUT TAPE 4,1.G,B,A,T 

1 F0RMATI3H G=F 1 0. 5 , 3X . 2HD=F 1 0 .5, 3X ,2 HA=F 10. 5. 3X, 2HT=F 1 0. 5/ ) 

WRITE OUTPUT T APE 4, 2 , T ,TEMP, SUM , DI F , J , SU! J ) 

2 FORMAT !3H T=F 10.5, 3X, 5HTEMP=E20 .6,3 X ,4HSUM=E20. 8, 3X,4HD I F=F 1 0.9 , 
13X,3HSU! I3,2H)=E20.8////) 

T=T+TD 

IF!T-TM)204,204,101 
GO TO 101 
10 STOP 10 
END 



*TEMP2 is the same except statement 300. 
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Slomax Subroutine (cont'd) 



<3 



SLOV=Slope 
T7=T(4)-T(2) 







SLOT 

=SL04 

r(1)=T(l+) 







T7 




•DIFP = 




DIF] 


F =■ 


SL04 


Cv 


SL04 


^L 05 




-S] 


L05 




SL03=SL02 

T(3)=T(2) 

SL02=SL04 

T(2)=T(1+) 

■SLQfaSLQlf 

(Outpu'O/ 
OUT 1 







SL01=SL02 
T(1)=T(2) 
SL02=SL04 
T(2)=T(4) 
-SLQ5= p LQ^ + . 

(Output^ 
OUT 1: 



C gf ) 
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SL03 

=SLOJ+ 

T(3)=0X4) 



(Output) 
OUT 27 

— pz — 
1 
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SUBROUTINE SLOMAX 

. DIKENSION LT ITLE ( 10) ,RR ( 126 ) »RI ( 123 ) ,T (10) , SLO( 10) 

COMMON RR.RI , T , SLO , G , B , A, P I , SLOPE ,M T , NT 1 , NT2, EP 1 til ,T1 

COMMON NT3.TM 

IF(SL0PE)4l,42,43 

41 WRITE OUTPUT TAPE 4t44 

44 F0RMAT(20H SLOPE 2 IS NEGATIVE/) 

GO TO 43 

42 WRITE OUTPUT TAPE 4,45 

45 FORMAT(26H CALCULATED SLOPE IS ZERO./) 

43 IF(NT1)5,6,6 
5 IF(NT)2,3,4 

2 K=1 
T(K)=T1 
SL01=SL0P£ 

, SLO(K)=SLOPE 
IF (NT3)39,38,39 

38 Tl=1 .*T(2)/2. 

GO TO 40 

39 Tl=(B+Tl)/2. 

40 K=K+1 
NT = 0 
RETURN 

3 T(K)=T1 
SL02=SL0PE 
SLC(K)=SLOPE 
SL05=SL02 

IF (SL02-SL01 ) 19,21,21 
19 T1=1 .*T( 1 )/4. 

NT=-1 
NT3 = 0 
PRINT 35 

35 F0RMAT(7H OUT 19/) 

RETURN 

21 T1=TM 
K=K+ 1 

, NT=1 
■ RETURN 

4 T(K)=T1 
SL03=SL0PE 

IF (SL02-SL03)23,22,22 
23 T1=4.*T(3) /3. 

NT=1 

PRINT 36 

36 F0RMAT(7H OUT 23/) 

RETURN 

22 SLO(K)=SLOPE 
K = K+1 

1 C1=SL01*(T(2)*T(2)-T(3)*T(3) ) 

C2=-SL02*(T( 1)«T(1)-T(3)*T(3)) 

C3=SL03*(T( 1 )*T( 1 )-T(2)«T(2) ) 

C4=SL01«(T(2)-T(3) ) 

C5=-SL02*(T( 1)-T(3)) 

C6 = SL03»(T( 1 )-T(2) ) 

T1=.5*(C1+C2+C3)/(C4+C5+C6) 

NT1 =0 

PRINT 30, K 
30 F0RMATI3H K=I3/) 
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C. 



RE TURN 

6 ER1=0.C 
SLCU=SLOPE 
T(K) = n 
SLO(K)=SLOPE 
T7=T(U)-T(2) 
IF(T7)7,8,9 

7 0IFF=SL0U-SLC5 
IF(DIFF)25,20,20 

25 SL01=SL04 
SLO( 1 )=SL04 
T( 1 )=T(4) 

PRINT il 

31 F0RRAT(7H OUT 25/) 

GO TO 1 

20 IF(DIFF-EPl) 14,14, 15 

8 CONTINUE 
NT2=1 
RETURN 

9 01 FF=SL04-SLC5 
IF (DIFF)27,26,26 

27 SLC3=SLC4 
SL0(3)=SLO4 
T(3)=T(4) 

PRINT 32 

32 F0RMAT(7H OUT 27/) 

GO TO 1 

26 IF(DIFF-EPl) 11,11,12 

11 CONTINUE 
NT2 = 1 
RETURN 

12 CONTINUE 
SL01=SL02 
SLC( 1 )=SL02 
T(1)=T(2) 

SL02=SL04 

SL0(2)=SL04 

T(2)=T(4) 

SL05=SL04 
PRINT 33 

33 F0RPAT(7H OUT 12/) 

GO TO 1 

14 CONTINUE 
NT2 = 1 
RETURN 

15 CONTINUE 
SL03=SL02 
SL0(3)=SL02 
T(3)=T(2) 

SL02=SL04 . 

SL0(2)=SL0U 
T(2 )=T( 4) 

SLC5=SL04 
PRINT 34 

34 F0RMAT(7H OUT 15/) 

GO TO 1 

END 
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SUBROUTINE AVER (S, A. II , J 1 , C I F , SUM ) 

DIMENSION SI200) , A( 12, 1C) 

12=11+1 
N1=J1-I 1-1 
DO 1 1=1,12 

A(I,1)=0.5*(S(N1)+S(N1+1) ) 

N1 = N1 + 1 . 

13=1 1-1 
N1 = I2 

DO 2 1 = 1 , 1 3 
N1=N1-1 
DO 2 J=1 ,N1 

A(J,I+l)=0.5»(A(J+l,I)+A(J,in 

DIF = ABSF( ( AINl , 13+1 )-A( N 1-1 , 13 + 1 ))/ AINl , 13+1 n 

SUM=A(N1 ,13+1 ) 

END 

SUBROUTINE GAUSSC ( 0, E , SUM ) 

DIMENSION XSI (24), WSI ( 24 ) , X ( 9 ) , Wl 9) , RR (1 28 ) , R I ( 1 28 ) ,T ( 1 0 ) ,SLO( 10) 
COMMON RR.RI ,T,SLO,G,B, A,P I,SL0PE,MT,NT1 ,NT2,EP1 , 1 1 ,T1 
COMMON NT3,TM 

IF(W)5,1,5 , 

W=1 .0 

XS I ( 16)=0. 1834346425 
XSI ( 171=0.5255324099 
XSI ( 18)=0. 7966664774 
XSI ( 19)=0. 9602098565 
XSI(20)=0.0 
XSI (21 )=0. 3242534234 
XSI (22)=0. 61 33714327 
XSI (23)=0. 836031 1073 
XSI (24)=0. 9681602395 
WSI ( 16)=0. 3626837834 
WSI ( 17)=0. 3137066459 
WSI ( 18)=0. 2223810345 
WS I ( 19)=0. 1012285363 
WSI (20)=0. 3302393550 
WSI (21 )=0. 3123470770 
WSI (22)=0. 2606106964 
WS I (23)=0. 1806481607 
WSI (24)=0. 0812743884 
Nl=8 
12 = 19 
I3=Nl/2 
K2 = I2 

A1=(E+D)*0.5 

A2=(E-D)*0.5 

I5=N1 

DO 35 1=1,13 

X( I )=A1-A2*XSI (K2) 

W( I )=A2*WSI (K2) 

X( I5)=A1 +A2+XSI (K2) 

W( I5)=W( I ) 

K2=K2-1 
15=15-1 
SUM=0.0 
DO 37 1=1, N1 
X1=X( I ) 

CALL EVAL35(X1 ,VR,VI ) 

SUM=SUM+VI*W( I ) 

RETURN 

END 
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SUBROUTINE SCOTS2 ( R. X , Y ) 

DIMENSION RM29) ,X( 129) tRR( 128) .RK 128) 
COMMON G,B,A,T,TM,TO,EPfLl,LO, I U mI ,RR, 
R( 1 ) = 1.0 
R(2)=B 

R(3)=-(B/A)*(G+1.) 

R(U)=-( (B**2)/A)*G 

X( 1 )=0.0 

X(2)*0.0 

X(3)=-(6/A)*Y 

X(4)=-( (8**2)/A)*Y 

NO=3 

MM=0 

IT=0 

CALL TOOTS2(R.X,NO,IT,MM,-.5,+.5) 

IF (RR{ 1 ) )4,4,S 
RRR=0. 

RII=0. , 

RRR=RR(1) ) 

RII=RI{1) 

RR(1)=RR(2) 

RI ( 1 )=RI (2) 

RR(2)=RR(3) 

RI (2)=RI (3) 

RR(3)=RRR 
RI(3)=RII 
GO TO 3 
IF (RR(2) )3,3f6 
RRR=RR( 2) 

RII=RI(2) 

RR(2)=RR(3) 

RI {2)=RI(3) 

RR(3)=RRR 
RI (3)=RII 
RETURN 
END 





, IT( 10) 
RI,PI 
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For use with SL02 and SL04. 



SUBROUTINE E VAL2S { RR , R I , G , B , Y , C 1 , , V I ) 

DIMENSION Z( 3) ,SR{ 4,4 ) ,SI ( 4,4 ) , E( 3) , S( 3) ,C( 3) tRR (128) ,R I M28) ,W(3) 
DO 100 1=1,3 

Z( I ) = ( (RR{ I ) »RR( I )-RI ( I )*R I ( I ) )/0)<-RR( I ) 

100 W( I ) = (2.*RR( I )*RI( I )/B) +RK I ) 

DOlOl 1=1,2 

D0102 J=2,3 

SR( I, J) = (Z( I )*Z( J) )-(W( I)*W( J) ) 

102 SI { I , J) = (Z{ I )*W( J) ) + (Z( J)*W( I ) ) 

101 CONTINUE 
E31=EXPF(RR(3)+RR( 1 ) ) 

S31=SINF(Rl (3)+RI( 1 ) ) ‘ 

C31=COSF{RI(3)+RI( 1 ) ) 

AR1=£31*C31 

AI1=£31*S31 

E3=EXPF{RR{3)) 

S3=SINF(RI (3)) 

C3=COSF(RI (3) ) 

DR1=E3*C3 
DI1=E5»S3 
E1=EXPF(RR{ 1 ) ) 

S1=SINF(RI { 1 n 

Cl l=COSF(RI { 1 ) ) 

ER1=E1*C11 

£I1=E1«S1 

SRN=SR(2,3)-SR{ 1 ,2 ) 

SIN=SI (2,3)-SI ( 1,2) 

SRDWSR(2,3)-SR{ 1,3) 

SID1=SI(2,3)-SI(1,3) 

SRD2=SR{ 1,3)-SR(1,2) 

SID2 = SI ( 1,5)-SI( 1,2) 

ANR=(SRN*AR1 )-(SIN*AI 1 ) 

ANI = (SIN*AR1 ) + {SRM»AIl ) 

DNR1=SRD1*DR1-SI01 *DI 1 
DNI 1=SID1»DR1+SRD1«0I 1 
0NR2=SRD2*ER1-SID2*£I 1 
0NI2=SID2*£R1+SR02»E11 
DNR=DNR1+DNR2 

DNI=0NI 1+DNI2 I 

ANUR=ANR*ONR+ANI*ONI 
ANUI=ANI*DNR-ANR*DNI 
DENOM=ONR*ONR+ONI.*DNI 
. VR=C1*{ ANUR/DENOM) 

VI=C1*{ ANUI/DENOM) 

END 
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W( I ) =(2. »RR( I ) *RI n ) /B) +Rl ( I ) 



102 

101 



103 



bbiol 1=1 f 2 

?R( i?j)=?zn )*z(j) )-(w( i)*w( j) ) 

si(i:j)=1z(i)*w(j))+(Z(J)*w(I)) 
CONTINUE ,,, nni^,, 
631=EXPF(RR(3) +RR( 1 ) ) 
S31=SINF(RI{3)+RI 1 
C31=C0SF(RI ( 3)+RI 1 ) ) 

E32=EXPF(RR(3)+RRI2) ) 
S32=SINF(RI(3)+RI(2 
C32=C0SF(RI(3 +RI 2 
E21=EXPF(RR(2)+RR 1)) 
C21=C0SF(RI (2)+RI 1 

S21=SINF(RI(2)+RIl 1) ) 

DO 103 1=1 »3 
E ( I )=eXPF( RR( I ) ) 

S( 1 )=SINF( RI ( I ) ) 

C( I )=COSF( RI { I ) ) 

AR1=E31*C31 
AI 1=E31*S31 
AR2=E21*C21 
AI2=E21*S21 
8Rl=AR2 



For use with SL02 and SL04 



when NTU less than 5 



BI 1=AI2 
^ BR2 = E32»C32 
BI2=E32*S32 
CR1=BR2 
Cl 1=BI2 
CR2=ARI 
CI2=AI 1 
AR12=AR1-AR2 
AI 12=AI 1-AI2 

A1R=(SR(2,3)»AR12)-(SI(2 
AI I=(SI (2,3)»AR12)+(SR( 2 
DR1=E(3)*C(3) 

DI 1=E(3)*S(3) 
DR2=E(2)*C(2) 
DI2=E(2)*S(2) 
DR12=DR1-0R2 
0112=011-012 
ER1=E(1 )*C( 1 ) ' 

EI1=E( 1 )*S( 1 ) 

ER2=DR1 
EI2=DI 1 
ER12=ER1-ER2 
El 12 = EI 1-EI2 
FR1=DR2 
FI 1=012 
FR2=ER1 
FI2=EI 1 
FR12=FR1-FR2 
FI 12=F1 1-FI2 
RR12=BR1-BR2 
BI 12=BI 1-BI2 
CR12=CR1-CR2 
Cl 12=CI 1-CI2 

B2R=(SR( 1,31*BR12)-(SI( 1 
B2I = (SI ( 1,3) «BR12) + ( SRI 1 
X3R=(SR(1,2)*CR12)-(SI( 1 
C3I=(SI ( 1,2 )*CR12)+(SR( 1 
D1R=(SR(2,3)*DR12)-(SI(2 
01 I = (SI (2,3)*0R12)+(SR(2 
E2R=(SR( 1,3)*ER12)-(SI I 1 
E2I=(S1 I 1,3)*ER12)+(SR( 1 
F3R=(SR( 1, 2)*FRl2)-( SII 1 
F3I=(SI I 1,2)*FR12)+(SR( 1 
ANR=A1R+B2R+C3R 
ANI=A1 I+B2I+C3I 
0NR=01R+E2R+F3R' 

ONI=01 I+E2I+F3I 

ANUR=ANR*ONR+ANI*ONI 

ANUR=ANR»ONR+ANI*ONI 

ANUI=ANI*ONR-ANR*ONI 

ANUI=ANI«ONR-ANR«ONI 

0EN0M=0NR»0NR+0NI*0NI 

0EN0M=0NR*0NR+0NI*0N1 

VR=C1*( ANUR/OENOM) 

VI=C1*( ANUI/DENOM) 

END 



,3)*AI 1 2) 
,3)*AI 12) 



,3)«BI 12) 
,3)»BI12) 
,2)*CI 12) 
,2)*CI12) 
,3)*0I12) 
,3)*0I12) 
,3)*EI12) 
,3)*EI12) 
,2)*FI12) 
,2)*FI12) 
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For use with TEMP2 atuf TEMP4 






■ SUBROUTINE E VAL 53S ( X 1 , VR , V I ) 

DIMENSION Z(3)»SR(4,4),SI(4,U),E(3),S(3),C{3)»RR(128),Rin28),W(3) 
1 ,R (129) ,X( 129) 

COMMON G,B,A,T,TK,TO,EP,Ll,LD,Il,N1,RRTRIfPI 

Y= ( PI/T) *X1 

CALL SOOTS2(R,X,Y) 

Cl=1.0 

DO 100 1=1 , 3 

Z( I ) = ( ( RR( I ) »RR( I )-RI ( I )»RI ( I ) ) /B ) + RR( I ) 

100 W( I ) = (2.»RRI'I)«RI( I )/B)+RI ( I ) 

D0101 1=1,2- 

D0102 J=2,3 

SR(I,J)=(2(I)*Z(J))-(W(I)*h(J)) 

102 SI ( I , J )=(Z( I )*W( J ) )+( Z( J ) »W( I ) ) 

101 CONTINUE 
r31=EXPF(RR(3)+RR( 1 ) ) 

S31=SINF(RI(3)+RI( 1) ) 

C31=C0SF(RI (3)+RI( 1 ) ) 

AR1=£3UC31. 

AI1=E31*S31 

E3 = EXPF(RR(3) ) 

S3 = SINF(RI( 3) ) 

C3 = C0SF(RI (3) ) 

DR1=E3*C3 
0I1=E3«-S3 
E1=£XPF(RR( 1 ) ) 

S1=SINF(RI ( 1 ) ) 

Cl 1=C0SF(RI( 1 ) ) 

ER1=E1«C11 

EI1=E1*S1 

SRN = SR(2,3)-SR( 1,2) 

SIN = SI(2,3)-SI( 1,2) 

SRD1=SR(2,3)-SR(1,3) 

SID1=SI (2,3)-SI(l,3) 

SRU2=SR( 1 ,3)-SR (1,2) 

SI02=SI( 1 ,5 )-SI (1,2) 

ANR=(SRN*AR1 )-( SIN*AI 1 ) 

ANI = (SIN<^AR1 ) + (SRN^*AIl ) 

0NR1=SR01 «DR1-S IC1*0I 1 
■ DNn=SI0UDR1+SRDl*0Il 
DNR2 = SR02*ER1-S in2»EI 1 
0Nl2 = SID2»ERl+SRD2«^En 
0NR=DNR1+DNR2 
DNI=DNI 1 +DNI2 
DENR=(G*DNR) -Y«-0NI 
0ENI=(Y*DNR) +G*DNI 
0ENCM= (OENR*DENR) + DEN UDENI 

VR=C1*C0SF( PI«X1 )* (ANR*CENR+ANI *DENI ) /OENOM 
VI =-Cl*SINF ( PI»X1 ) *( ANI »CENR-ANR*DENI )/DENOM 
END 
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For use with SL03‘ and SL05, R = 0. 



SUBROUTINE EVAL35 ( X 1 . VR , V I ) 

DIMENSION VR(20) ,VI( 20) ,SU(200),AV( 1 2 , 1 0 ) , R { 1 29 ) , X ( 1 29) ,ER(20C), 
1RR( 128) ,RI ( 128) ,T( 10) ,SLO( 10) 

COKKON RR. RI »T,SLO,G» B, A, P I , SLOPE t T , NT 1 , NT2 , EP 1 , 1 1 , T 1 
COMPON NT3,TM 
Y=(PI/T1)*X1 
Cl = l . 

01 = ( (G*G)+G+(Y»Y) )/( (G+1. )«(G+1 . )+( Y»Y) ) 

E1 = EXPF(-B*D1 ) 

F=Y/({G+1. )*(G+1.)+(Y*Y)) 

C=CCSF(B*F) 

S=SINF( B*F) 

VR = C1*E1*C*C0SF(PI*X1 ) 

VI*C1*E1*S*SINF( PI*X1 ) 

END 

CMO 



For use with TEMP2 and TEMP4, R = 0 



SUBROUTINE EVAL3T ( X 1 , VR . VI ) 

cT=i,o 

G1=G+1. 

DEN=G1*G1+Y*Y 
01 = (G*G+Y*Y ) /OEN 
F1=(G1»Y-G*Y)/DEN 
E1=EXPF(-B*D1 ) 

C*C0SF(B»F1) 

S=SINF( 8*F1 ) 

DENON=G*G+Y*Y 

VR=C1«C0SF (PI*X1 )*E1* (C*G+S*Y )/0EN0M 
V I =C1 *S INF ( P I*X 1 )*£!*( S*G-C«Y)/DENDM 
END 



S(3) ,C(3) ,RR( 128) ,R I ( 128) ,W(3) 
RRfRItPI 
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